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1  Introduction 


Two  directions  of  work  have  been  pursued  under  this  effort.  The  first  direc  ion  i  ^ 
with  further  extensions  of  our  geometric  techniques  for  spatial  and  surface  la}  lao  g 
applications  of  these  techniques  to  problems  in  high-frequency  electromagnetics.  le  sec¬ 
ond  direction  is  concerned  with  development  of  generalized  eikonal  and  transport  equations 
which  can  represent  multiphase  solutions  and  diffraction  effects.  The  first  approach  uses  di¬ 
rect  variational  and  geometrical  methods  for  building  numerically  high-frequency  solutions 
accounting  for  diffraction  effects.  It  is  applicable  to  geometrically  complex  objects  such 
as  aircraft,  ships,  satellites,  etc.  The  second  approach  is  based  on  solving  numerically  the 
generalized  eikonal  and  transport  equations  and  has  the  capabilities  to  capture  the  behavior 
of  solutions  near  caustic  regions  where  multiphase  solutions  may  develop. 


The  results  obtained  so  far  have  shown  that  both  directions  of  work  lead  to  interesting 
theoretical  results  and  numerical  codes  with  capabilities  for  solving  important  and  difficult 
high-frequency  wave  propagation  problems.  The  two  approaches  are  complementary  to  eacli 
other,  and  together  they  cover  a  large  variety  of  problems  in  electromagnetics. 

We  describe  now  the  results  in  more  detail. 


2  Calculation  of  geometric  quantities  required  for  com 
puting  currents  and  fields  due  to  creeping  waves 


The  work  on  this  task  involved  five  parts: 

1.  Investigation  of  efficient  means  for  estimating  the  required  length  of  continuation  of 
geodesics  into  shadow  region 

•>.  Development  of  fast  methods  for  determining  shadow  boundaries  suitable  for  UTl.) 
calculations 

3.  Testing  the  accuracy  of  the  algorithm  for  identifying  the  shadow  boundary  on  aircraft 
models 

4.  Upgrading  and  testing  of  the  algorithm  and  code  for  calculating  the  radii  of  curvature 
on  faceted  surfaces 

0.  Upgrading  and  testing  of  the  algorithm  and  code  for  building  initial  approximations 
to  geodesics. 
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In  order  to  calculate  currents  in  the  shadow  region  one  needs  to  construct  the  geodesics 
emanating  from  the  shadow  boundary  into  the  shadow  region.  Our  methods  for  computing 
geodesics  on  faceted  surfaces  are  very  fast.  In  particular,  the  number  of  operations  requiiec 
to  calculate  a  single  geodesic  of  length  not  exceeding  the  intrinsic  diameter  of  the  surlace  is 
at  most  of  order  0(N),  where  N  is  the  number  of  facets.  Usually,  however,  this  estimate  is 
overly  pessimistic  and  the  actual  calculation  of  a  geodesic  requires  processing  of  only  a  small 
fraction  of  of  facets,  therefore  reducing  significantly  the  computational  time. 

In  certain  cases,  important  for  electromagnetic  calculations,  one  needs  to  determine 
either  a  portion  or  an  entire  wavefront  of  a  creeping  wave  produced  by  a  source  on  the 
surface  or  in  space.  In  such  cases  the  wavefront  can  be  a  complicated  set  of  points  on 
the  surface.  An  accurate  determination  of  this  wavefront  may  require  tracing  a  dense  set 
of  geodesics  emanating  from  the  source  and  the  issue  of  computational  time  may  become 
critical.  Typically,  creeping  wavefronts  must  be  computed  in  situations  when  one  needs  to 
calculate  currents  produced  by  incident  waves  penetrating  into  the  shado\\  legion.  In  these 
cases  a  way  of  reducing  the  computational  time  needed  to  determine  the  required  geodesics 
is  to  construct  only  the  parts  of  geodesics  on  which  the  fields  are  above  a  certain  significant 
level.  This  is  based  on  the  fact  that,  according  to  the  results  by  V.  A.  Fock  and  other 
authors,  the  field  intensity  decays  as  the  wave  propagates  into  the  shadow  region.  Therefore, 
it  is  desirable  to  have  a  criteria  for  terminating  a  geodesic  once  the  field  intensity  is  below 
some  critical  value.  Application  of  such  a  criteria  will  allow  to  save  on  the  computational 
time  required  for  construction  of  geodesics. 

We  in\-estigated  this  approach  by  using  the  numerical  values  of  the  Fock  integral  as 
the  required  criteria.  W’e  have  developed  an  algorithm  and  code  for  calculating  the  Fock 
integral  along  an  arbitrary  surface  geodesic.  The  implemented  version  requires  that  the 
surlace  file  be  represented  as  a  collection  of  flat  triangular  facets.  However,  the  algorithm  is 
quite  general  and  can  be  implemented  with  other  surface  representations,  in  particular,  for 
surfaces  represented  by  Nonuniform  Rational  Bi-cubic  Splines  (NURBS).  Computationally, 
the  evaluation  of  the  Fock  integral  is  fast  and  the  results  obtained  so  far  show  that  the 
proposed  way  of  dealing  with  this  issue  is  satisfactory. 

.Another  necessary  step  in  computing  currents  and  fields  due  to  creeping  waves  is  the 
determination  of  the  shadow  boundary.  The  shadow  boundary  has  to  be  identified  in  two 
cases:  (a)  the  source  (or  the  field  point)  is  at  infinity  and  (b)  the  source  (or  the  field  point) 
is  at  a  finite  distance  from  the  surface.  For  the  case  (a)  we  developed  a  very  fast  (of  order 
log  A  .  where  A'  is  the  number  of  facets  in  a  surface)  algorithm  that  determines  the  shadow 
boundary  (boundaries)  on  an  arbitrary  flat-faceted  surface.  The  algorithm  is  based  on  a 
■■s])reading  spot"  strategy.  In  this  approach  we  pick  first  an  arbitrary  facet  (call  it  F)  and 
chc'ck  the  visibility  condition  for  this  facet.  If  the  facet  F  is  “visible”  to  the  source,  the 
program  checks  the  neighboring  facets  for  visibility  and,  if  visible,  adds  them  to  the  list  of 
visible  facets.  Once  all  the  neighboring  facets  of  F  are  examined,  the  facet  F  is  marked  and 
not  examined  any  more.  Then  all  the.  neighboring  facets  of  F  are  processed  in  the  same 
fashion.  This  process  continues  until  all  facets  that  are  visible  to  the  source  and  can  be 
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linked  with  F  through  a  chain  of  visible  facets  are  found. 

If  the  facet  F  is  invisible  to  the  source,  we  build  in  the  same  way  as  before  a  collection 
of  facets  not  visible  to  the  source.  As  a  result  of  this  procedure  the  entire  surface  is  divided 
into  the  lit  and  shadow  regions.  The  boundaries  of  this  regions  are  the  required  shadow 

boundaries. 

This  algorithm  has  been  coded  and  tested  on  various  surfaces,  including  aircraft 
models  such  as  VFY218,  Global  Hawk.  KC135,  and  F15.  The  code  showed  good  performance 
and  for  canonical  surface,  such  as  spheres,  ellipsoids,  paraboloids,  etc.,  the  computed  lesiills 
were  in  agreement  with  results  that  can  be  predicted  by  analytic  techniques. 

It  should  be  noted  that  at  this  stage  this  algorithm  does  not  identify  the  occluded 
parts  of  a  surface  and  therefore  a  part  of  the  surface  that  is  “blocked  from  the  souice  b_\ 
another  part  of  the  surface  will  still  be  identified  as  visible.  How^ever.  in  many  practical 
electromagnetic  problems  the  results  that  can  be  computed  with  the  already  de^•eloped 
algorithm  are  sufficient. 

In  the  case  (b)  when  the  source  (or  the  field  point)  is  at  a  finite  distance  from  the  sui- 
face  only  a  minor  modification  of  the  algorithm  is  required  to  identify  the  shadow  boundaries. 
This  modification  and  the  coding  have  been  implemented  and  tested. 


The  next  step  consists  in  constructing  a  sufficiently  dense  set  of  surface  geodesics 
tha1  originate  on  the  shadow  boundary  and  penetrate  into  the  shadow  region.  The  distances 
by  which  these  geodesics  are  continued  into  the  shadow  region  are  determined  according 
to  the  criteria  based  on  the  values  of  the  Fock  integral  as  described  earlier.  The  actual 
const  l  uction  of  the  geodesics  is  done  with  the  use  of  an  algorithm  that  we  developed  earlier. 
The  onlv  modification  that  is  needed  here  is  to  organize  the  calculations  of  geodesics  in  a 
way  suitable  for  dealing  with  many  “source”  points  located  on  the  shadow  boundary.  The 
issue  of  positioning  these  equivalent  sources  on  the  shadow  boundary  is  delicate  because  the 
shadow  boundary  may  (and  usually  will)  have  a  complicated  geometry.  We  are  working 
currently  on  this  task.  We  tested  the  algorithm  for  identifying  the  shadow  boundary  as 
rec(uired  for  high  frequency  calculations.  The  results  are  satisfactory  but  more  work  is 
needed  in  order  to  improve  the  accuracy.  The  difficulties  are  connected  with  the  fact  that, 
th('  t\'])ical  aircraft  geometry  is  quite  complex  and  the  usual  acceleration  strategies  needed 
for  fast  computation  of  the  shadow  boundary  rely  on  certain  trade-offs  that  impact  the 
accuracy.  We  are  continuing  working  on  this  issue. 

Se\eral  improvements  have  been  made  to  our  original  code  for  computing  radii  of 
curvature  of  surface  curves  on  faceted  surfaces.  These  improvements  enhance  the  accuracy  of 
the  calculations.  We  have  also  implemented  a  version  of  the  Dijkstra  algorithm  for  com])uting 
initial  approximations  to  surface  geodesics.  The  Dijkstra  algorithm  is  a  very  efficient  way 
for  computing  paths  on  faceted  surfaces  that  pass  through  surface  vertices.  In  our  scheme. 
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the  output  of  the  Dijkstra  algorithm  provides  initial  approximations  that  are  fed  into  oui 
optimization  code  that  determines  the  actual  geodesics.  Previously,  m  order  to  ^he 
required  initial  approximations  we  used  a  specially  developed  projection  algorithm.  While 
the  Dijkstra  algorithm  is  more  robust  than  our  projection  algorithm,  it  does  not  always 
capture  all  of  the  required  geodesics.  Consequently,  we  are  using  now  a  combined  scheme 
in  which  both,  our  original  projection  algorithm  and  the  Dijkstra  algorithm  are  employed. 
The  results  are  quite  satisfactory  in  terms  of  robustness  and  accuracy. 


3  Calculation  of  high  frequency  asymptotic  expan¬ 
sions 

\\b  have  used  the  FDTD  technique  applied  to  some  simple  geometries  in  order  to  validate 
the  methods  developed  in  this  project.  The  codes  for  the  FDTD  simulation  have  earlier  been 
developed  in  tw'o  and  three  space  dimensions.  They  have  now  been  adjusted  to  be  useful  tor 
validation. 

VVe  have  concentrated  on  the  development  of  generalized  eikonal  equations.  Two  types 
of  g('neralizations  have  been  studied.  One  is  the  derivation  of  partial  differential  equations 
which  can  represent  multiphase  solutions.  These  solutions  correspond  to  crossing  rays.  The 
otlK'f  is  the  inclusion  of  diffraction  effects. 

The  multiphase  technique  has  now  been  extended  to  the  time  periodic  case.  This 
allows  for  the  possibility  of  representing  more  crossing  rays  and  gives  a  computationally 
faslei-  algorithm.  This  method  does  hot  represent  diffraction  phenomena  but  gives  a  way  to 
coni])ute  the  ray  paths. 

The  equations  which  include  diffraction  effects  are  of  Schroedinger  equation  type. 
Thp\-  giv('  \-ery  good  representation  of  diffraction  phenomena  and  crossing  rays  as  long  as 
th('  angles  Iretween  the  rays  are  small.  We  have  developed  a  robust  numerical  method. 

We  obtained  a  new  set  of  partial  differential  equations  (PDFs)  which  can  be  seen 
as  a  generalization  of  the  classical  eikonal  and  transport  equations,  to  allow  for  solutions 
with  multiple  irhases.  The  traditional  geometrical  optics  pair  of  equations  suffer  from  the 
fact  that  the  class  of  physical  relevant  solutions  is  limited.  In  particular,  it  does  not  include 
solutions  with  multiple  phases,  corresponding  to  crossing  waves.  Our  objective  has  been 
to  generalize  these  e<iuations  to  accommodate  solutions  containing  more  than  one  phase. 
The  new  equations  are  based  on  the  same  high  frequency  approximation  of  the  scalar  wave 
equation  as  the  eikonal  and  the  transport  equations.  However,  they  also  incorporate  a  finite 
su])crposition  princijrle.  The  maximum  allowed  number  of  intersecting  waves  in  the  solution 
can  be  chosen  arbitrarily.  Irut  a  higher  number  means  that  a  larger  system  of  PDFs  must  be 
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solved.  The  PDEs  form  a  hyperbolic  system  of  conservation  la,ws  with  source  terms. 

Although  the  equations  are  only  weakly  hyperbolic,  and  thus  not  well  posed 
strong  sense,  several  examples  show  the  viability  of  solving  the  equations  numerically.  The 
technique  we  use  to  capture  multivalued  solutions  is  based  on  a  closure  assumption  tor  a 
system  of  equations  representing  the  moments. 

In  this  contract  we  thus  developed  a  middle  way  between  geometrical  optics  and  the 
kinetic  model.  It  is  a  high-frequency  approximation  through  which  the  whole  field  can  be 
solved.  ■  Moreover,  the  superposition  principle  holds  up  to  a  point;  the  maximum  allowed 
number  of  intersecting  waves  can  be  chosen  arbitrarily,  but  a  higher  number  means  that  a 
larger  system  of  PDEs  must  be  solved. 

The  starting  point  for  this  approach  is  the  transport  equation.  Instead  ol  solving 
the  full  equation  in  phase  space,  we  observe  that  when  /  is  of  a  simple  form  in  p.  we  can 
transform  (1)  to  a  finite  system  of  moment  equations  in  the  reduced  space  (xj).  analogously 
to  the  classical  approach  of  the  hydrodynamic  limit  from  a  kinetic  formulation. 

ft  +  v-VJ  +  cV^n-VJ  =  0  (1) 

In  particular,  we  are  interested  in  cases  where,  for  given  ;r  and  the  density  function  j  is 
nonzero  only  for  a  finite  number  of  p.  This  corresponds  to  a  finite  number  ol  rays  in  different 
directions  at  each  point. 

The  moment  ecjuations  are  derived  from  the  kinetic  model  for  high-frecpiency  wa\  es. 
They  are  equivalent  to  the  equations  of  geometrical  optics.  We  also  explored  some  theoretical 
issues  and  find  that  the  resulting  hyperbolic  equations  are  not  well-posed  in  the  strong  sense. 
Existence  of  solutions  of  unbounded  variation  is  indicated.  Next,  we  solved  these  equations 
for  one  and  two  phases.  The  standard  Lax-Friedrichs  method  gives  satisfactory  results.  Most 
elaborate,  and  less  viscous,  methods  like  the  Godunov  method  and  the  second-order  TV  I) 
Nessyahu-Tadmor  scheme,  although  converging  well  in  Li,  suffer  from  problems  locally  and 
converge  poorly  in  .  For  the  two-phase  system,  the  sensitivity  of  the  equations  is  more 
pronounced  and  consequently  it  is  harder  to  find  stable  numerical  methods.  .4fter  proper 
initialization,  the  ecjuations  can  however  be  solved  with  the  Lax-Friedrichs  method. 

The  research  on  multiphase  eikonal  equations  has  continued  by  testing  new  closure 
assumptions. 

The  original  closure  assumption  for  the  system  of  moment  equations  was  a  finite  sum 
of  delta  functions  in  the  space  of  ray  directions.  That  means  we  assumed  not  more  than  a 
finite  number  of  rays  at  each  point  in  space.  The  new  assumption  is  instead  a  finite  sum  of 
characteristic  functions  for  intervals.  The  endpoints  for  the  intervals  correspond  to  rays  in 
the  geometrical  optics  approximation. 


The  advantage  of  the  new  assumption  is  that  the  final  system  of  conservation  laAM> 
only  describes  the  ray  pattern  and  not  the  amplitudes.  It  appears  that  this  system  will  bo 
more  robust  at  caustics.  A  working  code  for  this  has  been  tested  for  3  superimposed  phases 
and  is  available  upon  request. 

In  addition,  the  FDTD  research  has  continued  in  two  directions.  The  parallelization 
of  the  code  has  been  improved  and  tested  on  SMP  nodes.  The  bandwidth  to  memory  gi\es 
the  limitations  in  the  computations.  A  version  of  the  code  is  available  upon  request. 

In  order  to  be  able  to  simulate  the  electromagnetic  field  near  small  geometric  details 
a  new  techniciue  for  wavelet  compression  of  the  difference  operators  has  been  tested  lor 
the  Helmholtz  equation.  Before  the  compression  the  difference  approximation  describes  the 
geometry  in  detail  on  a  refined  grid.  The  compression  results  in  a  coarser  grid  method  which 
still  approximates  the  fine  geometric  details. 

.  This  is  part  of  work  involving  a  systematic  technique  for  the  derivation  of  subgrid  scale 
models  in  the  numerical  solution  of  partial  differential  equations.  The  technique  is  based  on 
Haar  wavelet  projections  of  the  discrete  operator  followed  by  a  sparse  approximation.  .A.S 
numerical  testing  suggests,  the  resulting  numerical  method  will  accurately  represent  subgrid 
scale  phenomena  on  a  coarse  grid.  Applications  to  wave  propagation  in  materials  with 
siibgi  id  inhomogeneities  are  demonstrated  below. 


We  consider  the  Helmholtz  equation. 


for  </(.r)  >  0. 


-L 


dxi 


-  k^u{x)  =  0. 


\V('  discretize  (2)  in  one  dimension,  with  the  boundary  conditions  w(0)  =  1  and 

(/'(!)=  0. 

—  -^A+aiA-iii  —  k^ui  0.  r  =  l,...,n. 

uo  =  2td0)  —  ui,  (3) 

We  use  k  =  2-  and  (i(x)  =  1 1  if  0.45  <  r  <  0.551  otherwise  and  we  take  n  =  256  and  use 
three  lioinogenizat ions.  We  get 

Lfii  =  (Lj  -  k^iyu^O.  (4) 

Truncation  is  performed  on  I,  (or  H j)  and  not  on  Lj.  The  result  is  shown  in  Figure  1. 
\V('  see  that  Helmholtz  equation  gives  results  similar  to  those  of  the  model  equation.  Band 
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projection  is  more  efficient  than  truncation  and  working  on  H,  the  subgrid  model,  is  more 

efficient  than  working  on  I,,  the  Schur  complement. 

We  consider  the  two-dimensional  version  of  (2)  with  periodic  boundary  conditions  in 
the  y-direction.  and  u(0,y)  =  1,  =  0  at  the  left  and  right  boundaries  respectiveh 

This  corresponds  to  a  plane  time-harmonic  wave  of  amplitude  one  entering  the  computational 
domain  from  the  left  and  flowing  out  at  the  right.  The  discretization  that  we  use  is 


—AlauAiuu  -  k^Uic 


iJ^l. 


=  UiA. 


This  leads  to  the  matrix  equation 

=  U.FeVj+u  n  =  m2^-^\  (b) 

where  m  is  a  positive  integer  and  Tj+i  is  homogenized  following  our  theory  for  two-dimensional 
problems. 

.A.S  an  example  we  choose  the  aix.y)  shown  in  Figure  2.  which  models  a  wall  with  a^ 
small  slit  where  the  incoming  wave  can  pass  through.  W  ith  k  =  Sir  and  n  =  48.  we  obtained 
the  results  presented  in  Figure  3. 


The  structure  of  the  operator  after  one  homogenization  step  is  shown  in  Figure  4. 
Truncation  of  this  operator  according  to  (5)-(6)  gives  the  result  shown  in  Figure  .5  for  various 
values  of  The  case  ;/  =  9  corresponds  to  a  compression  to  approximately  7%  of  the  original 
size. 
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